#pragma rtGlobals=1		// Use modern global access method.

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------
Menu "Spectrometry"
	"Maximum Entropy", maxEntGUI()
end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function resFunction(xWidGuass, xWidLorent, imardp)
Variable xWidGuass, xWidLorent
Wave imardp
Variable p1,p2,suma, sumb

Make/O/N=(dimsize(imardp,0)) blardp, blarsqdp, guassTemp, lorenTemp, func

xWidGuass = xWidGuass/(2*sqrt(2*ln(2)))

guassTemp = exp(-0.5*((x-dimsize(imardp,0)/2)/(xWidGuass))^2)  
lorenTemp = (1/pi)*(xWidLorent/2)/((x-dimsize(imardp,0)/2)^2 + (xWidLorent/2)^2)

If(xWidGuass != 0 && xWidLorent != 0)
	Duplicate/O guassTemp, func
	Convolve/A lorenTemp, func
elseif(xWidGuass != 0 && xWidLorent == 0)
	func = guassTemp
elseif(xWidGuass == 0 && xWidLorent != 0)
	func = lorenTemp
endif

blardp = func

blarsqdp = blardp^2

suma = sum(blardp)
sumb = sum(blarsqdp)

blardp = blardp/suma
blarsqdp = blarsqdp/sumb

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function setupArrays(imardp, guess)
Wave imardp,guess

Make/O/N=(dimsize(imardp,0),dimsize(imardp,1)) maparr, errardp, vari, d2chi

maparr = guess

Variable nconst,p1,p2

nconst = sum(imardp)
Variable/G nconstG = nconst

Wavestats/Q maparr

If(V_min < 2*10^-9)
	setzero("maparr")
endif
 
normal("maparr", nconst)

Variable numberOfPnts
If(dimsize(maparr, 1) == 0)
	numberOfPnts = numpnts(maparr)
elseif(dimsize(maparr, 1) != 0)
	numberOfPnts = (dimsize(maparr,0)*dimsize(maparr,1))
endif

Variable meanVal = nconst/numberOfPnts

errardp[][] =  (maparr[p][q] - meanVal)^2/(numberOfPnts-1)
setzero("vari")


p1=0
p2=0

Do
	p2 = 0
	Do 
		If(errardp[p1][p2] < 5)
			errardp[p1][p2] = 4
		endif
		p2+=1
	while(p2 < dimsize(errardp,1))
	p1+=1
while(p1 < dimsize(errardp,0))

vari = 2/errardp
d2chi = vari
//dconv(vari,blarsqdp,"d2chi")

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function setzero(zarrayStr)
String zarrayStr
Wave zarray = $zarrayStr
Variable p1=0,p2=0

Do
	p2 = 0
	Do
		if(zarray[p1][p2] < 1*10^-5)
			zarray[p1][p2] = 1*10^-5
		endif
		p2+=1
	while(p2 < dimsize(zarray,1))
	p1+=1
while(p1 < dimsize(zarray,0))


end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function normal(narrayStr, nconst)
String narrayStr
Variable nconst

Wave narray = $narrayStr
Variable suma,p1,p2

Do
	p2 = 0
	Do
		suma = suma + narray[p1][p2]
		p2+=1
	while(p2 < dimsize(narray,1))
	p1+=1
while(p1 < dimsize(narray, 0))

narray = narray*nconst/suma

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function dconv(carray,psf,conresStr)
Wave carray,psf
String conresStr
Make/O/N=(dimsize(carray,0),dimsize(carray,1)) $conresStr
Make/O/N=(dimsize(carray,0),dimsize(carray,0)) funcMat
Wave conres = $conresStr, funcMat
Variable p1,p2,p3,p4, bb

conres = carray
funcMat[][] = psf[p + dimsize(psf,0)/2 - q]

MatrixOp/O conres = funcMat x carray

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function getchi(rawdata, consol, varianceArray)
Wave rawdata, consol, varianceArray
Variable chisq, n, redchi, p1, p2, s
Make/O/N=(dimsize(rawdata,0),dimsize(rawdata,1)) chiar

chiar = (((rawdata - consol)^2)/(varianceArray))
chisq = sum(chiar)
n = (numpnts(rawdata))
redchi = chisq/(sum(rawdata)- numpnts(rawdata)-1)

Make/O/N=3 returnWave = {chisq, n, redchi}

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function mainIterLoop(imardp)
Wave imardp
NVAR lambda = lambda,plag = plag
NVAR nconst = nconstG, oldchi = oldchi
Wave blardp, blarsqdp, d2chi, dupmaparr, errardp, inputMat, maparr, vari

Variable sdist = 0, cdist = 0, p1,p2
Make/O/N=(dimsize(maparr,0), dimsize(maparr,1)) workarr, entropyarr, gradientS, gradientC
workarr = maparr

entropyarr[][] = (maparr[p][q]-maparr[p-1][q])^2/(2*errardp[p][q])*exp(-(maparr[p][q]-maparr[p-1][q])^2/(2*errardp[p][q]))
gradientS[][] = (maparr[p][q]-maparr[p-1][q])/(errardp[p][q])*(exp(-(maparr[p][q]-maparr[p-1][q])^2/(2*errardp[p][q]))*((maparr[p][q]-maparr[p-1][q])^2/(2*errardp[p][q]))+1)

Do
	p2 = 0
	Do
		sdist = sdist + maparr[p1][p2]*((maparr[p1][p2]-maparr[p1-1][p2])/(errardp[p1][p2])*(exp(-(maparr[p1][p2]-maparr[p1-1][p2])^2/(2*errardp[p1][p2]))*((maparr[p1][p2]-maparr[p1-1][p2])^2/(2*errardp[p1][p2]))+1))^2
		p2+=1
	while(p2 < dimsize(maparr,1))
	p1+=1
while(p1 < dimsize(maparr,0))

sdist = sdist^(1/2)

dconv(workarr,blardp,"mapblur")

Wave mapblur

Make/O/N=(dimsize(maparr,0),dimsize(maparr,1)) dir1, dir2, dir3, condir1, condir2, condir3, rworkarr

p1 = 0
p2 = 0

Do
	p2 = 0
	Do
		dir1[p1][p2] = -1*maparr[p1][p2]*((maparr[p1][p2]-maparr[p1-1][p2])/errardp[p1][p2]*(exp(-(maparr[p1][p2]-maparr[p1-1][p2])^2/(2*errardp[p1][p2]))*((maparr[p1][p2]-maparr[p1-1][p2])^2/(2*errardp[p1][p2]))+1))
		p2+=1
	while(p2 < dimsize(maparr,1))
	p1+=1
while(p1 < dimsize(maparr,0))

MatrixOp/O dir2 = 2*(imardp - mapblur)/errardp

rworkarr = dir2

dconv(rworkarr,blardp,"dir2")

MatrixOp/O sumz = sum(maparr*(dir2*dir2))
cdist = sumz[0]
MatrixOp/O dir2 = -maparr*dir2

cdist = cdist^(1/2)
gradientC = 2*(imardp - mapblur)/errardp


MatrixOP/O dir3 = (d2chi*dir1)/sdist - (d2chi*dir2)/cdist)

Variable s1,s2,s3,s11,s22,s33,s12,s13,s23,d11,d22,d33,d13,d12,d23

MatrixOp/O sumz = sum(dir1 - entropyarr*gradientS)
s1 = sumz[0]
MatrixOp/O sumz = sum(dir2 -entropyarr*gradientC)
s2 = sumz[0]
MatrixOp/O sumz = sum(dir3 -((entropyarr*gradientS*vari)/sdist)-(entropyarr*gradientC*vari)/cdist))
s3 = sumz[0]
MatrixOp/O sumz = sum((dir1*dir1)/(2*maparr))
s11 = sumz[0]
MatrixOp/O sumz = sum((dir2*dir2)/(2*maparr))
s22 = sumz[0]
MatrixOp/O sumz = sum((dir3*dir3)/(2*maparr))
s33 = sumz[0]
MatrixOp/O sumz = sum((dir1*dir2)/(maparr))
s12 = sumz[0]
MatrixOp/O sumz = sum((dir1*dir3)/(maparr))
s13 = sumz[0]
MatrixOp/O sumz = sum((dir2*dir3)/(maparr))
s23 = sumz[0]
MatrixOp/O sumz = sum((dir1*dir1)/2)
d11 = sumz[0]
MatrixOp/O sumz = sum((dir2*dir2)/2)
d22 = sumz[0]
MatrixOp/O sumz = sum((dir3*dir3)/2)
d33 = sumz[0]
MatrixOp/O sumz = sum((dir1*dir2))
d12 = sumz[0]
MatrixOp/O sumz = sum((dir1*dir3))
d13 = sumz[0]
MatrixOp/O sumz = sum((dir2*dir3))
d23 = sumz[0]

Make/O/N=(dimsize(dir3,0),dimsize(dir3,1)) tworkarr

workarr = dir1/maparr
rworkarr = dir2/maparr
tworkarr = dir3/maparr

dconv(workarr,blardp,"condir1")
dconv(rworkarr,blardp,"condir2")
dconv(tworkarr,blardp,"condir3")


Variable c1, c2, c3, c11, c22, c33, c12, c13, c23

MatrixOp/O sumz = sum((imardp - mapblur)/errardp*condir1)
c1 = sumz[0]
MatrixOp/O sumz = sum((imardp - mapblur)/errardp*condir2)
c2 = sumz[0]
MatrixOp/O sumz = sum((imardp - mapblur)/errardp*condir3)
c3 = sumz[0]
MatrixOp/O sumz = sum((condir1*condir1)/errardp)
c11 = sumz[0]
MatrixOp/O sumz = sum((condir2*condir2)/errardp)
c22 = sumz[0]
MatrixOp/O sumz = sum((condir3*condir3)/errardp)
c33 = sumz[0]
MatrixOp/O sumz = sum((condir1*condir2)/errardp)
c12 = sumz[0]
MatrixOp/O sumz = sum((condir1*condir3)/errardp)
c13 = sumz[0]
MatrixOp/O sumz = sum((condir2*condir3)/errardp)
c23 = sumz[0]

c1 = -2*c1
c2 = -2*c2
c3 = -2*c3
c11 = .5*c11
c22 = .5*c22
c33 = .5*c33

Make/O/N=(3,3) matrix

Make/O/N=(3,1) conVect, solVect

matrix[0][0] = s11 - lambda*c11 - plag*d11
matrix[0][1] = s12 - lambda*c12 - plag*d12
matrix[0][2] = s13 - lambda*c13 - plag*d13
conVect[0][0] = lambda*c1 - s1
matrix[1][0] = s12 - lambda*c12 - plag*d12
matrix[1][1] = s22 - lambda*c22 - plag*d22
matrix[1][2] = s23 - lambda*c23 - plag*d23
conVect[1][0] = lambda*c2 - s2
matrix[2][0] = s13 - lambda*c13 - plag*d13
matrix[2][1] = s23 - lambda*c23 - plag*d23
matrix[2][2] = s33 - lambda*c33 - plag*d33
conVect[2][0] = lambda*c3 - s3

MatrixOp/O solVect = inv(matrix) x conVect

Variable alpha1, alpha2, alpha3, change, entropy

alpha1 = solVect[0]
alpha2 = solVect[1]
alpha3 = solVect[2]

p1 = 0
p2 = 0

NoteBook Maximum_Entropy text="Test value = " + num2str(maparr[50][50]) + "\r"

MatrixOP/O changeMat = alpha1*dir1 + alpha2*dir2 + alpha3*dir3

Do
	p2 = 0
	Do
		if (abs(changeMat[p1][p2]/maparr[p1][p2]) > 1)
			changeMat[p1][p2] = changeMat[p1][p2]/abs(changeMat[p1][p2])*maparr[p1][p2]*.99
		endif
		entropy = entropy + (maparr[p1][p2]-maparr[p1-1][p2])^2/(2*errardp[p1][p2])*exp(-(maparr[p1][p2]-maparr[p1-1][p2])^2/(2*errardp[p1][p2]))
		p2+=1
	while(p2 < dimsize(maparr,1))
	p1+=1
while(p1 < dimsize(maparr,0))

maparr = maparr + changeMat
workarr = maparr

NoteBook Maximum_Entropy text="Test value now = " + num2str(maparr[50][50]) + "\r"

setZero("maparr")

normal("workarr", nconst)
normal("maparr", nconst)

dconv(workarr,blardp,"mapblur")

normal("mapblur", nconst)

getchi(imardp, mapblur, errardp)

Wave returnWave

Variable chisq = returnWave[0], npts = returnWave[1], redchi = returnWave[2]

NoteBook Maximum_Entropy text="Chi-squared = " + num2str(chisq) + " on " + num2str(npts) + " pixels.\r"
NoteBook Maximum_Entropy text="Reduced Chi-squared = " + num2str(redchi) + "\r"
NoteBook Maximum_Entropy text="Entropy = " + num2str(entropy) + "\r"
NoteBook Maximum_Entropy text="Lambda = " + num2str(lambda) + "\r"
NoteBook Maximum_Entropy text="Plag = " + num2str(plag) + "\r\r"

Make/O/N=(dimsize(maparr,0), dimsize(maparr,1)) oldmap

if(redchi > oldchi)
	maparr = oldmap
	plag = plag*10
elseif(redchi <= oldchi)
	oldmap = maparr
	oldchi = redchi
endif

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function maximumEntropy(inputWaveStr, guessWaveStr, energyWaveStr, xWidGuassEV, xWidLorentEV, lambdaIn, plagIn, iterations, smoothing)
String inputWaveStr, guessWaveStr, energyWaveStr
Variable xWidGuassEV,xWidLorentEV,lambdaIn, plagIn, iterations, smoothing
Variable/G oldchi = 1*10^9, lambda = lambdaIn, plag = plagIn
Variable xWidGuass, xWidLorent
Variable p1
Wave inputWave = $inputWaveStr, guessWave = $guessWaveStr, energyWave = $energyWaveStr
SVAR saveDF = root:Packages:saveDF

WaveStats/Q inputWave
If(V_max < 700 && dimsize(inputWave, 1) == 0)
	inputWave = (inputWave - V_min)/(V_max - V_min)
	inputWave = inputWave*750
endif

WaveStats/Q energyWave
Variable delta = (V_max - V_min)/(numpnts(energyWave))
xWidGuass = xWidGuassEV/delta
xWidLorent = xWidLorentEV/delta

resFunction(xWidGuass,xWidLorent,$inputWaveStr)
setupArrays(inputWave,guessWave)

if(WinType("Maximum_Entropy") != 0)
	DoWindow/K  Maximum_Entropy
	DoWindow/K  MaxEnt_Graph
endif
NewNoteBook/N=Maximum_Entropy/F=0/W=(10,10,190,260) as "Maximum_Entropy"

	Do
		NoteBook Maximum_Entropy text="Iteration Number: " + num2str((p1+1)) +"\r"
		mainIterLoop(inputWave)
		If(dimsize(inputWave,1) != 0)
			convertTo1D("maparr")
			convertTo1D(inputWaveStr)
		endif
		If(p1 == 0)
			If(dimsize(inputWave,1) != 0)
				String nameHolder = inputWaveStr + "_1D"
				Wave inputWave_1D = $nameHolder
				Display/N=MaxEnt_Graph/W=(200,100,800,500)  inputWave_1D,maparr_1D vs energyWave
				ModifyGraph rgb(maparr_1D)=(0,0,0)
			else
				Display/N=MaxEnt_Graph/W=(200,100,800,500) inputWave,maparr vs energyWave
				ModifyGraph rgb(maparr)=(0,0,0)
			endif
		endif
		DoUpdate
		If(plag >= 10^6)
			break
		endif
		
	if(dimsize(inputWave, 1) == 0)
//		Smooth smoothing, maparr
	endif
			

		p1+=1
	while(p1 < iterations)
	
If(dimsize(inputWave,1) != 0)	
	convertTo1D(inputWaveStr)
	Duplicate/O maparr_1D,maparr_1D_smth
	Smooth smoothing, maparr_1D_smth
elseif(dimsize(inputWave, 1) == 0)
	Duplicate/O maparr,maparr_smth;DelayUpdate
	Smooth smoothing, maparr_smth
endif
//curvatureCorrect("maparr","curvature")
//curvatureCorrect(inputWaveStr,"curvature")
DoWindow/K MaxEnt_Graph
If(dimsize(inputWave,1) != 0)
	String nameHolder1 = inputWaveStr + "_1D"
	Wave inputWave_1D = $nameHolder1
	Display/N=MaxEnt_Graph/W=(200,100,800,500)  inputWave_1D,maparr_1D_smth vs energyWave
	ModifyGraph rgb(maparr_1D_smth)=(0,0,0)
else
	Display/N=MaxEnt_Graph/W=(200,100,800,500) inputWave,maparr_smth vs energyWave
	ModifyGraph rgb(maparr_smth)=(0,0,0)
endif
SetDataFolder saveDF
end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------
Function maxEntGUI()
NewDataFolder/O root:Packages
Variable/G root:Packages:xWidGuassEV = 0,root:Packages:xWidLorentEV = 0,root:Packages:lambdaIn = 1000, root:Packages:plagIn = 10, root:Packages:iterations = 100, root:Packages:smoothing = 2
NVAR xWidGuassEV = root:Packages:xWidGuassEV, xWidLorentEV = root:Packages:xWidLorentEV, lambdaIn = root:Packages:lambdaIn, plagIn = root:Packages:plagIn, iterations = root:Packages:iterations, smoothing = root:Packages:smoothing
String/G root:Packages:inputWaveStr ="", root:Packages:guessWaveStr ="", root:Packages:energyWaveStr=""
SVAR  inputWaveStr = root:Packages:inputWaveStr, guessWaveStr = root:Packages:guessWaveStr, energyWaveStr = root:Packages:energyWaveStr
String/G root:Packages:waveSelectionReturnList = ""
DoWindow/K/Z maximumEntropyWindow
	NewPanel /W=(150,77,584,306) /N=maximumEntropyWindow as "Maximum Entropy"
	ModifyPanel frameStyle=3
	Button inputWaveButton,pos={10,45},size={50,20},title="Browse", proc=browseButton_Proc
	Button guessWaveButton,pos={10,105},size={50,20},title="Browse", proc=browseButton_Proc
	Button energyWaveButton,pos={10,165},size={50,20},title="Browse", proc=browseButton_Proc
	SetVariable guassianWidth,pos={274,15},size={145,16},title="Gaussian Width   ", value=xWidGuassEV
	SetVariable guassianWidth,limits={-inf,inf,0}
	SetVariable lorentzianWidth,pos={274,45},size={145,16},title="Lorentzian Width ", value=xWidLorentEV
	SetVariable lorentzianWidth,limits={-inf,inf,0}
	SetVariable lambdaValue,pos={274,75},size={145,16},title="Lambda               ", value=lambdaIn
	SetVariable lambdaValue,limits={-inf,inf,0}
	SetVariable plagValue,pos={274,105},size={145,16},title="Plag                     ",value=plagIn
	SetVariable plagValue,limits={-inf,inf,0}
	SetVariable iterationsValue,pos={274,135},size={145,16},title="Iterations             ",value=iterations
	SetVariable iterationsValue,limits={-inf,inf,0}
	SetVariable smoothingValue,pos={274,165},size={145,16},title="Smoothing           ",value=smoothing
	SetVariable smoothingValue,limits={-inf,inf,0}
	SetVariable inputWaveValue,pos={10,15},size={250,16},title="Input Wave   ",value=inputWaveStr
	SetVariable inputWaveValue,limits={-inf,inf,0},noedit=1
	SetVariable guessWaveValue,pos={10,75},size={250,16},title="Guess Wave ",value=guessWaveStr
	SetVariable guessWaveValue,limits={-inf,inf,0},noedit= 1
	SetVariable energyWaveValue,pos={10,135},size={250,16},title="Energy Wave",value=energyWaveStr
	SetVariable energyWaveValue,limits={-inf,inf,0},noedit= 1
	Button maxEntButton,pos={310,195},size={50,20},title="Execute", proc=executeButton_Proc
	Button cancelButton,pos={375,195},size={50,20},title="Cancel", proc=genericCancel_Proc
	ModifyPanel fixedSize=1
	SetWindow maximumEntropyWindow, hook(windowCancel_Hook) = windowHookFunction_Cancel
end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------
Function executeButton_Proc(ctrlName) : ButtonControl
String ctrlName
NVAR xWidGuassEV = root:Packages:xWidGuassEV, xWidLorentEV = root:Packages:xWidLorentEV, lambdaIn = root:Packages:lambdaIn, plagIn = root:Packages:plagIn, iterations = root:Packages:iterations, smoothing = root:Packages:smoothing
SVAR  inputWaveStr = root:Packages:inputWaveStr, guessWaveStr = root:Packages:guessWaveStr, energyWaveStr = root:Packages:energyWaveStr
SVAR  waveSelectionReturnList = root:Packages:waveSelectionReturnList
String/G root:Packages:pathtoWaveFolder = ParseFilePath(1, Stringfromlist(0,waveSelectionReturnList), ":", 1, 0)
SVAR pathtoWaveFolder = root:Packages:pathtoWaveFolder
String/G root:Packages:saveDF
SVAR saveDF = root:Packages:saveDF
saveDF = GetDataFolder(1)
SetDataFolder $pathtoWaveFolder
maximumEntropy(inputWaveStr, guessWaveStr, energyWaveStr, xWidGuassEV, xWidLorentEV, lambdaIn, plagIn, iterations, smoothing)
end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function browseButton_Proc(ctrlName) :ButtonControl
String ctrlName
SVAR  inputWaveStr = root:Packages:inputWaveStr, guessWaveStr = root:Packages:guessWaveStr, energyWaveStr = root:Packages:energyWaveStr
SVAR  waveSelectionReturnList = root:Packages:waveSelectionReturnList
If(stringMatch(ctrlName,"inputWaveButton"))
	genericWaveSelection("Select a Wave", 100, 375, 100, 300, 0, 0, 200, 250)
	PauseForUser Wave_Browser
	inputWaveStr = ParseFilePath(0, Stringfromlist(0,waveSelectionReturnList), ":", 1, 0)
elseif(stringMatch(ctrlName,"guessWaveButton"))
	genericWaveSelection("Select a Wave", 100, 375, 100, 300, 0, 0, 200, 250)
	PauseForUser Wave_Browser
	guessWaveStr = ParseFilePath(0, Stringfromlist(0,waveSelectionReturnList), ":", 1, 0)
elseif(stringMatch(ctrlName,"energyWaveButton"))
	genericWaveSelection("Select a Wave", 100, 375, 100, 300, 0, 0, 200, 250)
	PauseForUser Wave_Browser
	energyWaveStr = ParseFilePath(0, Stringfromlist(0,waveSelectionReturnList), ":", 1, 0)
endif
end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------
Function convertTo1D(inputWaveStr)
String inputWaveStr
Wave w1 = $inputWaveStr
Variable suma, p1, p2

String nameHolder

nameHolder = inputWaveStr + "_1D"

Make/O/N=(dimsize(w1,0)) $nameHolder 
Wave w2 = $nameHolder

Do
	p2 = 0
	suma = 0
	Do
		suma = suma + w1[p1][p2]
		p2+=1
	while(p2 < dimsize(w1, 1))
	w2[p1] = suma
	p1+=1
while(p1 < dimsize(w1, 0))

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------
Function curvatureCorrect(inputWaveStr, curvatureWaveStr)
String inputWaveStr, curvatureWaveStr
Wave inputWave = $inputWaveStr, curvatureWave = $curvatureWaveStr
String nameHolder = inputWaveStr + "_1D"
Make/O/N=(dimsize(inputWave, 0)) $nameHolder
Wave w2 = $nameHolder, Channel

//splitMat(inputWaveStr)

Variable p1

Do
	nameHolder = inputWaveStr + num2str(p1+1)
	Wave w1 = $nameHolder
	MatrixOp/O w3 = col(inputWave, p1)
	Variable p2 = 0
	Do
		inputWave[p2][p1] = interp(Channel[p2] - curvatureWave[p1], Channel, w3)
		p2+=1
	while(p2 < dimsize(inputWave,0))
	p1+=1
while(p1 < dimsize(inputWave, 1))

convertTo1D(inputWaveStr)

KillWaves w3

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function entropyDerivatives(inputWaveStr)
String inputWaveStr
Wave inputWave = $inputWaveStr, vari

Make/O/N=(dimsize(inputWave,0),dimsize(inputWave,1)) entJacobian, entropyFunc
Make/O/N=(dimsize(inputWave,0),dimsize(inputWave,0),dimsize(inputWave,1)) entHessian

Variable p1,p2,p3

Do
	p2=0
	Do
		entropyFunc[p1][p2] = (inputWave[p1][p2]-inputWave[p1-1][p2])^2/(2*vari[p1][p2])*exp(-(inputWave[p1][p2]-inputWave[p1-1][p2])^2/(2*vari[p1][p2]))
		entJacobian[p1][p2] = (((inputWave[p1][p2]-inputWave[p1-1][p2])/vari[p1][p2])*(exp(-((inputWave[p1][p2]-inputWave[p1-1][p2])^2)/(2*vari[p1][p2]))*((inputWave[p1][p2]-inputWave[p1-1][p2])^2/(2*vari[p1][p2]))+1))
		p2+=1
	while(p2 < dimsize(inputWave,1))
	p1+=1
while(p1 < dimsize(inputWave,0))

p1 = 0
p2 = 0
p3 = 0

Do
	p2 = 0
	Do
		p3 = 0
		Do
			If(p1 != p2 && p1-1 != p2)
				entHessian[p1][p2][p3] = 0
			elseif(p1 == p2)
				entHessian[p1][p2][p3] = ((inputWave[p1][p3] - inputWave[p1-1][p3])/(vari[p1][p3]))*(exp(-(inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3]))*((inputWave[p1][p3]-inputWave[p1-1][p3])/vari[p1][p3])+((inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3])+1)*(-(inputWave[p1][p3]-inputWave[p1-1][p3])/(vari[p1][p3]))*exp(-(inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3])))
				entHessian[p1][p2][p3] = entHessian[p1][p2][p3] + (1/vari[p1][p3])*(exp(-(inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3]))*((inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3])+1))
			elseif(p1-1 == p2)
				entHessian[p1][p2][p3] = -((inputWave[p1][p3] - inputWave[p1-1][p3])/(vari[p1][p3]))*(exp(-(inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3]))*((inputWave[p1][p3]-inputWave[p1-1][p3])/vari[p1][p3])+((inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3])+1)*(-(inputWave[p1][p3]-inputWave[p1-1][p3])/(vari[p1][p3]))*exp(-(inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3])))
				entHessian[p1][p2][p3] = entHessian[p1][p2][p3] - (1/vari[p1][p3])*(exp(-(inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3]))*((inputWave[p1][p3]-inputWave[p1-1][p3])^2/(2*vari[p1][p3])+1))
			endif
			p3+=1
		while(p3 < dimsize(entHessian,2))
		p2+=1
	while(p2 < dimsize(entHessian,1))
	p1+=1
while(p1 < dimsize(entHessian,0))

end

//---------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------

Function makeGuessWave(inputWaveStr)
String inputWaveStr
Wave inputWave = $inputWaveStr
KillWaves/Z guessWave
Make/O/N=(dimsize(inputWave,0)) guessWave

Variable p1, V_flag, V_PeakLoc
WaveStats/Q inputWave
Do
	If(V_PeakLoc != dimsize(inputWave, 0)-1)
		FindPeak/B=0/R=[V_PeakLoc+1]/M=(.2*V_max)/P/Q inputWave
		guessWave[V_PeakLoc] = V_PeakVal
	endif
while(V_PeakLoc <= dimsize(inputWave, 0)-1)

end